R Notebook

R Notebook covering data generated from healthy samples for Figure 1, Supplementary figure 1 and 2.

LOAD ALL PACKAGES

library(tidyverse)

Registered S3 methods overwritten by 'dbplyr':
  method         from
  print.tbl_lazy     
  print.tbl_sql      
Registered S3 method overwritten by 'cli':
  method     from         
  print.boxx spatstat.geom
── Attaching packages ─────────────────────────────────────────────────────────────────────────────── tidyverse 1.3.1 ──
✓ ggplot2 3.3.5     ✓ purrr   0.3.4
✓ tibble  3.1.6     ✓ dplyr   1.0.7
✓ tidyr   1.1.4     ✓ stringr 1.4.0
✓ readr   2.1.0     ✓ forcats 0.5.1
── Conflicts ────────────────────────────────────────────────────────────────────────────────── tidyverse_conflicts() ──
x dplyr::filter() masks stats::filter()
x dplyr::lag()    masks stats::lag()

library(Seurat)

Attaching SeuratObject

library(cowplot)
library(RColorBrewer)
library(pheatmap)
library(scales)

Attaching package: ‘scales’

The following object is masked from ‘package:purrr’:

    discard

The following object is masked from ‘package:readr’:

    col_factor

Import functions from the custom R script

source("SPATIAL_FUNCTIONS.R")

package ‘clusterProfiler’ was built under R version 4.1.1
clusterProfiler v4.0.5  For help: https://yulab-smu.top/biomedical-knowledge-mining-book/

If you use clusterProfiler in published research, please cite:
T Wu, E Hu, S Xu, M Chen, P Guo, Z Dai, T Feng, L Zhou, W Tang, L Zhan, X Fu, S Liu, X Bo, and G Yu. clusterProfiler 4.0: A universal enrichment tool for interpreting omics data. The Innovation. 2021, 2(3):100141. doi: 10.1016/j.xinn.2021.100141

Attaching package: ‘clusterProfiler’

The following object is masked from ‘package:purrr’:

    simplify

The following object is masked from ‘package:stats’:

    filter

Loading required package: AnnotationDbi
Loading required package: stats4
Loading required package: BiocGenerics
Loading required package: parallel

Attaching package: ‘BiocGenerics’

The following objects are masked from ‘package:parallel’:

    clusterApply, clusterApplyLB, clusterCall, clusterEvalQ, clusterExport, clusterMap, parApply,
    parCapply, parLapply, parLapplyLB, parRapply, parSapply, parSapplyLB

The following objects are masked from ‘package:dplyr’:

    combine, intersect, setdiff, union

The following objects are masked from ‘package:stats’:

    IQR, mad, sd, var, xtabs

The following objects are masked from ‘package:base’:

    anyDuplicated, append, as.data.frame, basename, cbind, colnames, dirname, do.call, duplicated, eval,
    evalq, Filter, Find, get, grep, grepl, intersect, is.unsorted, lapply, Map, mapply, match, mget, order,
    paste, pmax, pmax.int, pmin, pmin.int, Position, rank, rbind, Reduce, rownames, sapply, setdiff, sort,
    table, tapply, union, unique, unsplit, which.max, which.min

Loading required package: Biobase
Welcome to Bioconductor

    Vignettes contain introductory material; view with 'browseVignettes()'. To cite Bioconductor, see
    'citation("Biobase")', and for packages 'citation("pkgname")'.

Loading required package: IRanges
Loading required package: S4Vectors
package ‘S4Vectors’ was built under R version 4.1.1
Attaching package: ‘S4Vectors’

The following object is masked from ‘package:clusterProfiler’:

    rename

The following objects are masked from ‘package:dplyr’:

    first, rename

The following object is masked from ‘package:tidyr’:

    expand

The following objects are masked from ‘package:base’:

    expand.grid, I, unname


Attaching package: ‘IRanges’

The following object is masked from ‘package:clusterProfiler’:

    slice

The following objects are masked from ‘package:dplyr’:

    collapse, desc, slice

The following object is masked from ‘package:purrr’:

    reduce


Attaching package: ‘AnnotationDbi’

The following object is masked from ‘package:clusterProfiler’:

    select

The following object is masked from ‘package:dplyr’:

    select

IMPORTING THE HEATLTHY SAMPLES

## HEALTHY SAMPLES
## HEALTHY MALE SKIN 1 & 2
HEALTHY.Male.s1 <- Load10X_Spatial(data.dir ="../../FOURTH_RUN/SAMPLE A1/",slice="ST-HM-1-R1")

'giveCsparse' has been deprecated; setting 'repr = "T"' for youInvalid name supplied, making object name syntactically valid. New object name is ST.HM.1.R1; see ?make.names for more details on syntax validityKeys should be one or more alphanumeric characters followed by an underscore, setting key from st.hm.1.r1_ to sthm1r1_

HEALTHY.Male.s2 <- Load10X_Spatial(data.dir ="../../SEVENTH RUN/ST-HM-1/",slice="ST-HM-1-R2")

'giveCsparse' has been deprecated; setting 'repr = "T"' for youInvalid name supplied, making object name syntactically valid. New object name is ST.HM.1.R2; see ?make.names for more details on syntax validityKeys should be one or more alphanumeric characters followed by an underscore, setting key from st.hm.1.r2_ to sthm1r2_

## HEALTHY FEMALE SKIN 1 & 2 (FROM THE SAME DONOR)
HEALTHY.Female1.s1 <- Load10X_Spatial(data.dir ="../../FOURTH_RUN/SAMPLE B1/",slice="ST-HF-1-R1")

'giveCsparse' has been deprecated; setting 'repr = "T"' for youInvalid name supplied, making object name syntactically valid. New object name is ST.HF.1.R1; see ?make.names for more details on syntax validityKeys should be one or more alphanumeric characters followed by an underscore, setting key from st.hf.1.r1_ to sthf1r1_

HEALTHY.Female1.s2 <- Load10X_Spatial(data.dir ="../../FIFTH_RUN/SAMPLE D1/",slice="ST-HF-1-R2")

'giveCsparse' has been deprecated; setting 'repr = "T"' for youInvalid name supplied, making object name syntactically valid. New object name is ST.HF.1.R2; see ?make.names for more details on syntax validityKeys should be one or more alphanumeric characters followed by an underscore, setting key from st.hf.1.r2_ to sthf1r2_

HEALTHY.Female1.s3 <- Load10X_Spatial(data.dir ="../../SEVENTH RUN/ST-HF-1/",slice="ST-HF-1-R3")

'giveCsparse' has been deprecated; setting 'repr = "T"' for youInvalid name supplied, making object name syntactically valid. New object name is ST.HF.1.R3; see ?make.names for more details on syntax validityKeys should be one or more alphanumeric characters followed by an underscore, setting key from st.hf.1.r3_ to sthf1r3_

## HEALTHY FEMALE
HEALTHY.Female2.s1 <- Load10X_Spatial(data.dir ="../../SEVENTH RUN/ST-HF-2E/",slice="ST-HF-2-R1")

'giveCsparse' has been deprecated; setting 'repr = "T"' for youInvalid name supplied, making object name syntactically valid. New object name is ST.HF.2.R1; see ?make.names for more details on syntax validityKeys should be one or more alphanumeric characters followed by an underscore, setting key from st.hf.2.r1_ to sthf2r1_

HEALTHY.Female2.s2 <- Load10X_Spatial(data.dir ="../../SEVENTH RUN/ST-HF-2F/",slice="ST-HF-2-R2")

'giveCsparse' has been deprecated; setting 'repr = "T"' for youInvalid name supplied, making object name syntactically valid. New object name is ST.HF.2.R2; see ?make.names for more details on syntax validityKeys should be one or more alphanumeric characters followed by an underscore, setting key from st.hf.2.r2_ to sthf2r2_

Import the spots to be removed from loupe browser.

# HEALTHY FEMALE SKIN 2
#REMOVE SPOTS

remove.spots <- read.csv(file="../../CLOUPE_FILES/FIFTH_RUN/D1/OUTLIERS.csv")
subset_spots <- Cells(HEALTHY.Female1.s2)[which((!(rownames(HEALTHY.Female1.s2@meta.data) %in% remove.spots$Barcode)))]
HEALTHY.Female1.s2.clean <- subset(HEALTHY.Female1.s2,cells=subset_spots)

#SAMPLE IDS

HEALTHY.Male.s1$sample.id<- "HV1_S1_R1" 
HEALTHY.Male.s2$sample.id<- "HV1_S1_R2"

HEALTHY.Female1.s1$sample.id <- "HV2_S1_R1"
HEALTHY.Female1.s2.clean$sample.id <- "HV2_S1_R2"
HEALTHY.Female1.s3$sample.id <- "HV2_S2"

HEALTHY.Female2.s1$sample.id <- "HV3_S1"
HEALTHY.Female2.s2$sample.id <- "HV3_S2"

SAMPLE IDS (version 2) - (Full name + Biopsy site)

TRUNK and FA (Forearm)

HEALTHY.Male.s1$orig.ident <- "Healthy Volunteer 1 - TRUNK R1" 
HEALTHY.Male.s2$orig.ident <-"Healthy Volunteer 1 - TRUNK R2"

HEALTHY.Female1.s1$orig.ident <- "Healthy Volunteer 2 - FA R1"
HEALTHY.Female1.s2.clean$orig.ident <- "Healthy Volunteer 2 - FA R2"
HEALTHY.Female1.s3$orig.ident <- "Healthy Volunteer 2 - TRUNK"

HEALTHY.Female2.s1$orig.ident <- "Healthy Volunteer 3 - FA R1"
HEALTHY.Female2.s2$orig.ident <- "Healthy Volunteer 3 - TRUNK"

Compile the samples into one list.

Healthy_Samples <- c(HEALTHY.Male.s1,HEALTHY.Male.s2,HEALTHY.Female1.s1,HEALTHY.Female1.s2.clean,HEALTHY.Female1.s3,HEALTHY.Female2.s1,HEALTHY.Female2.s2)

Generating spatial plots without QC filtering

for (x in Healthy_Samples){
  st_plot(x)
}

QC plots for healthy samples

for (x in Healthy_Samples){
  st_plot_QC(x)
}

QC scatter plots

for (x in Healthy_Samples){
  st_scatter_QC(x)
}

Filtering all samples to remove spots with less than 200 genes expressed

i <- 1
while(i <= length(Healthy_Samples)){
  filtered_data <- st_filter_by_genes(st.data = Healthy_Samples[[i]],x = 200)
  Healthy_Samples[[i]] <- filtered_data
  i <- i+1
}

Batch correction of all samples using Seurat Anchor integeration

new.skin.combined <- st_combine(Healthy_Samples,ndim = 20,res = 0.6)

ADDING ADDITIONAL META-DATA

meta.data_skin <- new.skin.combined@meta.data$orig.ident %>% as.data.frame()
write.csv(meta.data_skin,file="HEALTHY_SAMPLES_ST_META_DATA.csv")

IMPORTING ST META-DATA UPDATED VERSION

meta.data_skin_updated <- read.csv(file="HEALTHY_SAMPLES_ST_META_DATA.csv")

ADDING META-DATA to Seurat object

new.skin.combined@meta.data$orig.ident <- meta.data_skin_updated$orig.ident
new.skin.combined@meta.data$sample.id <- meta.data_skin_updated$sample.id

Determining marker genes for batch corrected healthy samples + heatmap

new.skin.combined.markers <- FindAllMarkers(new.skin.combined, only.pos = TRUE, min.pct = 0.25, logfc.threshold = 0.25,assay = "SCT")
DefaultAssay(new.skin.combined) <- "SCT"
top10 <- new.skin.combined.markers %>%
    group_by(cluster) %>%
    filter(p_val_adj<0.05) %>%
    top_n(n = 10, wt = avg_log2FC) %>%
    filter(gene %in% VariableFeatures(new.skin.combined_V2))
pdf(width = 20,height=15,file = "HEATMAP_HEALTHY_SAMPLES_CLUSTERS_ONLY.pdf")
print(DoHeatmap(new.skin.combined, features = top10$gene,assay = "SCT",group.colors = custom_colors) + NoLegend())
dev.off()

Assign cluster ids

new.cluster.ids <- c("0 Mixed","1 Dermis",
"2 Dermis Adipose","3 Epidermis 1","4 Dermis Connective Tissue","5 Upper follicle (UF) and Perifollicular dermis","6 Epidermis 2","7 Dermis 2","8 Dermis Smooth Muscle","9 Hair follicle","10 Dermis Lymphatics","11 Dermis vasculature")
names(new.cluster.ids) <- levels(new.skin.combined)
new.skin.combined <- RenameIdents(new.skin.combined, new.cluster.ids)
new.skin.combined <- StashIdent(new.skin.combined, save.name = "Spatial.regions")

Defining and viewing color scheme for UMAP and spatial plot

custom_colors <- c("#B4DFFC","#6EAB3D","#FFD700","#A020F0","#FFA500","#AEDD3C","#595959","#D2AF81FF","#3A49FC","#FF0000","#A86F3D","#A18FAA")
show_col(custom_colors)

General UMAP with defined color scheme

UMAP plot split by samples

DimPlot(new.skin.combined,split.by = "sample.id")

## SAVE THE RESULTS
save(new.skin.combined,file = "HEALTHY_SKIN_SAMPLES_ST.RData")

Spatial plot for HV2_S2_R1 /HEALTHY_FEMALE_SAMPLE_1

SAVE OUTPUT TO PDF

pdf(width = 8,height=10,file = "HEALTHY_FEMALE_SAMPLE_1_R2_SPATIAL_PLOT.pdf")
SpatialDimPlot(new.skin.combined,images="HealthyFemale.1.S2",cols = custom_colors,pt.size.factor = 2.4)
dev.off()

VIEW SPATIAL PLOT HERE

SpatialDimPlot(new.skin.combined,images="HealthyFemale.1.S2",cols = custom_colors,pt.size.factor = 2.4)

Generate Spatial plot for all samples

images <- Images(new.skin.combined) %>% as.vector()
for(x in images){
  pdf(width = 8,height=10,file = paste(x,"_SPATIAL_PLOT.pdf"))
  print(SpatialDimPlot(new.skin.combined,images=x,cols = custom_colors,pt.size.factor = 2))
  dev.off()
}

Color scheme for MIA heatmap

col.pal <- RColorBrewer::brewer.pal(9, "OrRd")

MIA - scRNA dataset (Travis et al)

SINGLE CELL DATA FROM TRAVIS PAPER

load("../../TRAVIS_scRNA_DATA/SkinSeuratTotal.Rdata")

# META DATA FROM TRAVIS 
meta.data <- read.csv("../../TRAVIS_scRNA_DATA/Cell_Level_Metadata.csv")

#AddMetaData(SkinSeuratTotal,metadata = meta.data$Specific_CellType,col.name = "Specific_Celltype")
skin.meta.data <- as.data.frame(SkinSeuratTotal@meta.data) %>% rownames_to_column("CellID")
inner_join(x=skin.meta.data,y=meta.data,by="CellID")

SkinSeuratTotal@meta.data <- inner_join(x=skin.meta.data,y=meta.data,by="CellID") %>% column_to_rownames("CellID")

SkinSeuratTotal@meta.data$Celltype <- meta.data$CellType
SkinSeuratTotal@meta.data$Condition <- meta.data$Condition
SkinSeuratTotal@meta.data$Celltype_specific <- meta.data$Specific
SkinSeuratTotal@meta.data$SampleID <- meta.data$SampleCondition

#SUBSET DATA TO ONLY HEALTHY SAMPLES
travis_NM_skin_scRNA<- subset(SkinSeuratTotal,subset = Condition == c("Normal"))

## PROCESS SC RNA DATA
travis_NM_skin_scRNA <- NormalizeData(travis_NM_skin_scRNA)

Performing log-normalization
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|

travis_NM_skin_scRNA <- FindVariableFeatures(travis_NM_skin_scRNA, selection.method = "vst", nfeatures = 2000)

Calculating gene variances
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Calculating feature variances of standardized and clipped values
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|

travis_NM_skin_scRNA <- ScaleData(travis_NM_skin_scRNA)

Centering and scaling data matrix

  |                                                                                                                      
  |                                                                                                                |   0%
  |                                                                                                                      
  |========================================================                                                        |  50%
  |                                                                                                                      
  |================================================================================================================| 100%

travis_NM_skin_scRNA <- RunPCA(travis_NM_skin_scRNA, features = VariableFeatures(object = travis_NM_skin_scRNA))

PC_ 1 
Positive:  SPARC, FTH1, COL6A2, DCN, FSTL1, PPAP2B, C1S, SERPINF1, COL6A1, COL1A2 
       C1R, COL6A3, IGFBP7, CTSK, CCDC80, TIMP2, VCAN, SPARCL1, FBLN1, MFAP4 
       AEBP1, FBLN2, SLIT3, IGFBP4, IGFBP5, ICAM1, MMP2, APOD, FN1, CFH 
Negative:  DMKN, SCEL, SBSN, KRTDAP, DSC1, CALML5, SPINK5, LYPD3, DSG1, TMEM45A 
       CSTA, KRT2, CASP14, TACSTD2, KLK7, IVL, SERPINA12, KRT10, CLDN4, POF1B 
       KRT80, ELOVL4, DHCR24, PERP, LY6G6C, GLTP, PKP1, MAL2, JUP, SLURP1 
PC_ 2 
Positive:  LCP1, LAPTM5, RGS1, SRGN, HLA-DQB1, TYROBP, HLA-DQA1, HLA-DPB1, GPR183, HLA-DRA 
       HLA-DRB1, HLA-DPA1, HLA-DRB5, CD83, FCER1A, CD74, HLA-DQA2, ARL4C, HLA-DQB2, HLA-DMA 
       PTPRC, GNA13, LST1, HLA-DMB, PLEK, CD86, TNFRSF1B, CST7, HLA-DOA, ITGB2 
Negative:  DCN, COL6A2, COL1A2, C1S, C1R, COL6A1, CTSK, CFD, COL6A3, FBLN1 
       FSTL1, MFAP4, TIMP2, CCDC80, VCAN, FBLN2, AEBP1, MMP2, IGFBP5, SLIT3 
       PPAP2B, DPT, SPARC, LUM, OLFML3, FN1, COL3A1, CFH, PCOLCE, APOD 
PC_ 3 
Positive:  CST3, HLA-DQA1, HLA-DPB1, HLA-DRA, HLA-DPA1, HLA-DRB1, HLA-DQB1, HLA-DRB5, C15orf48, CD74 
       HLA-DQA2, HLA-DQB2, LST1, HLA-DMB, HLA-DMA, CPVL, CD83, AIF1, PLEK, CTSS 
       TYROBP, CD86, CTSH, FCER1A, FAM49A, SERPINF1, TSPAN33, HLA-DOA, CSF1R, IFI30 
Negative:  KRT5, KRT14, TNS4, COL17A1, S100A2, KRT15, DSG3, AQP3, IGFBP3, POSTN 
       DSP, ITGA6, SFN, LAMB4, SYT8, BNC1, DST, LAMC2, ITGA2, FGFBP1 
       LAMA3, NRG1, PHLDB2, PERP, MET, EDN1, CDCP1, LAMB3, SERPINB2, IMPA2 
PC_ 4 
Positive:  PLVAP, CD93, ELTD1, DARC, TM4SF1, ECSCR, TSPAN7, AC011526.1, CLDN5, VWF 
       SELE, A2M, EMCN, S1PR1, C2CD4B, CDH5, IFI27, KDR, ERG, ESAM 
       GIMAP5, CSF3, FLT1, GNG11, PTPRB, SOX17, CNKSR3, TACR1, RCAN1, GIMAP6 
Negative:  ARL4C, DCN, PTPRC, KRT5, TNS4, DST, KRT14, LCP1, LAPTM5, CTSK 
       FBLN1, THBS2, CD69, CFD, COL1A2, SAMSN1, AQP3, MXRA5, SERPINF1, RUNX3 
       SYTL3, MFAP4, RGS1, C1S, POSTN, S100A2, DSG3, COL17A1, C1R, CDCP1 
PC_ 5 
Positive:  CPA3, TPSAB1, CTSG, IL1RL1, GATA2, KIT, HDC, GCSAML, HPGD, LAX1 
       SRGN, CPM, TPSD1, PAG1, BATF, SYTL3, CD69, CMA1, TNIK, SAMSN1 
       VWA5A, CALB2, RGS13, SDPR, HEY1, SLC2A3, CHN2, GPR65, HPGDS, TMOD1 
Negative:  KRT5, TNS4, KRT14, IL1R2, AQP3, DSG3, PLEK2, DSP, SFN, COL17A1 
       S100A2, SLC2A1, IGFBP3, HLA-DRA, HLA-DQA1, DST, CDCP1, CLDN1, HLA-DPA1, LAMB3 
       HLA-DPB1, HLA-DRB5, HLA-DQB2, FGFBP1, PERP, POSTN, HLA-DRB1, KRT15, CD74, HLA-DMA 

travis_NM_skin_scRNA <- FindNeighbors(travis_NM_skin_scRNA, dims = 1:40)

Computing nearest neighbor graph
Computing SNN

travis_NM_skin_scRNA <- FindClusters(travis_NM_skin_scRNA, resolution = 1)

Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck

Number of nodes: 4840
Number of edges: 159696

Running Louvain algorithm...

0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|

Maximum modularity in 10 random starts: 0.9054
Number of communities: 24
Elapsed time: 0 seconds

travis_NM_skin_scRNA <- RunUMAP(travis_NM_skin_scRNA, dims = 1:40)

The default method for RunUMAP has changed from calling Python UMAP via reticulate to the R-native UWOT using the cosine metric
To use Python UMAP via reticulate, set umap.method to 'umap-learn' and metric to 'correlation'
This message will be shown once per session12:58:23 UMAP embedding parameters a = 0.9922 b = 1.112
12:58:23 Read 4840 rows and found 40 numeric columns
12:58:23 Using Annoy for neighbor search, n_neighbors = 30
12:58:24 Building Annoy index with metric = cosine, n_trees = 50
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
12:58:24 Writing NN index file to temp file /var/folders/8p/2clwbcfn719f5l1kjsyx4v6s2r9kgj/T//RtmpdfWGXA/file1aa32e12e645
12:58:24 Searching Annoy index using 1 thread, search_k = 3000
12:58:25 Annoy recall = 100%
12:58:26 Commencing smooth kNN distance calibration using 1 thread
12:58:26 Initializing from normalized Laplacian + noise
12:58:27 Commencing optimization for 500 epochs, with 196850 positive edges
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
12:58:33 Optimization finished

Normal skin - scRNA data from Travis paper

Idents(travis_NM_skin_scRNA) <- "Specific_CellType"
DimPlot(travis_NM_skin_scRNA,pt.size = 2.5,label = TRUE)

Idents(travis_NM_skin_scRNA) <- "CellType"
DimPlot(travis_NM_skin_scRNA,pt.size = 2.5)

save(travis_NM_skin_scRNA,file = "../../TRAVIS_scRNA_DATA/travis_NM_skin_scRNA.RDS")

single_cell.markers <- FindAllMarkers(travis_NM_skin_scRNA,assay = "RNA",logfc.threshold = 0.25)

filtered_single_cell.markers <- single_cell.markers %>% filter(p_val_adj<=0.05) %>% group_by(cluster) %>% top_n(n =300,wt = avg_log2FC) %>% filter(avg_log2FC>0.25)
filtered_spatial_markers <- new.skin.combined.markers %>% filter(p_val_adj<=0.05) %>% group_by(cluster) %>% top_n(n =300,wt = avg_log2FC) %>% filter(avg_log2FC>0.25)

Figure - 2B

##INTERSECT GENES BETWEEN scRNA and Spatial data
st.genes <- unique(rownames(new.skin.combined@assays$Spatial@counts))
sc.genes <- unique(rownames(travis_NM_skin_scRNA@assays$RNA@counts))
all.genes.scrna_and_spt <- unique(intersect(sc.genes,st.genes))

MIA_results <- MIA(total_genes = length(all.genes.scrna_and_spt),single_cell.markers = filtered_single_cell.markers,spatial.markers = filtered_spatial_markers)

E.data <- MIA_results %>% column_to_rownames("cluster")
E.data <- E.data[,order(colnames(E.data))]
pheatmap(E.data,cluster_cols = FALSE,cluster_rows = FALSE,fontsize=15,color = col.pal)

immune_only.E.data <- E.data[,c("T cell-1","T cell-2","T cell-3","Myeloid-1","Langerhans","Mast-1","Mast-2")]
pdf(file = "MIA_regions_IMMUNE_CELLS.pdf",width = 10,height = 4)
pheatmap(immune_only.E.data,cluster_cols = FALSE,cluster_rows = FALSE,fontsize=15,color = col.pal)
dev.off()

structure_only.E_data <- E.data[,c("Fibroblast-1","Fibroblast-2","Fibroblast-4","Fibroblast-5","HairFollicle","Keratinocyte-2","Keratinocyte-3","Keratinocyte-5","Keratinocyte-6","Keratinocyte-8","Sebocyte","Venule-1","Venule-2","Venule-3","Venule-4","VSMC-1")]
pdf(file = "MIA_regions_STRUCTURAL_CELLS.pdf",width = 10,height = 4)
pheatmap(structure_only.E_data,cluster_cols = FALSE,cluster_rows = FALSE,fontsize=15,color = col.pal)
dev.off()

NOW CALCULATING ENRICHMENT USING NEW DATASET (Travis et al)

Seurat based cell type enrichment

# RUN THIS ON BIG PURPLE 
skin_reference.2 <- SCTransform(travis_travis_NM_skin_scRNA, ncells = 3000, verbose = FALSE) %>% RunPCA(verbose = FALSE) %>% 
    RunUMAP(dims = 1:40)
anchors <- FindTransferAnchors(reference = skin_reference.2, query = new.skin.combined, normalization.method = "SCT")
predictions.assay <- TransferData(anchorset = anchors, refdata = skin_reference.2$Specific_CellType, prediction.assay = TRUE, 
    weight.reduction = new.skin.combined[["pca"]],dims = 1:40)
new.skin.combined[["predictions_travis_data"]] <- predictions.assay

CELL TYPES FOR FIGURE 1 - Travis et al. single cell dataset (Only using the healthy subset of scRNA data)

Cell_types <- c("Fibroblast-1","HairFollicle","Keratinocyte-5","Keratinocyte-3","Sebocyte","VSMC-1","Melanocyte")
DefaultAssay(new.skin.combined) <- "predictions_travis_data"
for (x in Cell_types){
  pdf(file = paste("HEALTHY_FEMALE_1_R2_CELL_TYPE_",x,".pdf"),width = 10,height = 15)
  print(SpatialFeaturePlot(new.skin.combined, features = c(x), pt.size.factor = 2.5, ncol = 2, crop = TRUE,images = "HealthyFemale.1.S2") + scale_fill_gradientn(colors=cols,limits = c(0,1))) 
  dev.off()
}

PERCENTAGE PLOT FOR DIFFERENT SPATIAL REGIONS

Calculated on a per sample basis

Regions.df <- table(new.skin.combined@meta.data$Spatial.regions,new.skin.combined@meta.data$orig.ident) %>% as.data.frame() %>% dplyr::rename(Spatial_Region=Var1,Sample=Var2) 

black.bold.16.text <- element_text(face = "bold", color = "black", size = 14,angle = 90, vjust = 0.5, hjust=1)
brks <- c(0, 0.25, 0.5, 0.75, 1)

pdf(file="PERCENTAGE_COMPOSTION_PLOT_HEALTHY_SAMPLES.pdf",height = 14,width = 8)
ggplot(Regions.df,aes(x=Sample,y=Freq,fill=Spatial_Region)) + geom_bar(stat="identity",position="fill") + scale_fill_manual(values =custom_colors)  + ggplot2::theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
panel.background = element_blank(),axis.text.x =black.bold.16.text) + scale_y_continuous(breaks = brks, labels = scales::percent(brks)) + ylab("% Composition of Sample by Spatial Region / Clusters")
dev.off()

Cell counts for travis healthy subset data

cell.counts <- as.data.frame(table(travis_NM_skin_scRNA@meta.data$Specific_CellType))
write.csv(cell.counts,file="../../TRAVIS_scRNA_DATA/CELL_COUNTS_TRAVIS_HEALTHY_SKIN_DATA.csv")

Filtered out cell types with less than 50 cells

filtered_cell.counts <- cell.counts %>% filter(Freq>50) 
Idents(travis_NM_skin_scRNA) <- "Specific_CellType"
travis_healthy_data.subset <- subset(x = travis_NM_skin_scRNA, idents = c(filtered_cell.counts$Var1))

DimPlot(travis_healthy_data.subset,pt.size = 2.5)

Plotting final travis healthy data UMAP

pdf(width = 12,height=8,file = "UMAP_TRAVIS_DATA_HEALTHY_SAMPLES_ONLY.pdf")
DimPlot(travis_healthy_data.subset,pt.size = 3.5)
dev.off()

PATHWAY ANALYSIS FOR SPATIAL HEALTHY DATA - PER CLUSTER BASIS

for(i in seq(0,11)){
  cluster <- new.skin.combined.markers %>% filter(cluster==i)
  GO_PATHWAYS_ENRICH(cluster,ONT="BP",OUTPUT = i)
  GO_PATHWAYS_ENRICH(cluster,ONT="CC",OUTPUT = i)
  GO_PATHWAYS_ENRICH(cluster,ONT="MF",OUTPUT = i)
  KEGG_PATHWAYS(cluster,OUTPUT = i)
}

Number of genes used in MIA - (Needed for Venn diagram)

cluster.5 <- filtered_spatial_markers %>% filter(cluster==5)
Myeloid.1 <- filtered_single_cell.markers %>% filter(cluster=="Myeloid-1")

length(intersect(cluster.5$gene,Myeloid.1$gene))

cluster.6 <- filtered_spatial_markers %>% filter(cluster==6)
Keratinocyte.5 <- filtered_single_cell.markers %>% filter(cluster=="Keratinocyte-5")

length(intersect(cluster.6$gene,Keratinocyte.5$gene))

LOAD DATA FROM HANNIFA PAPER

load("../../scRNA_DATA/HEALTHY_SC_RNA_PROCESSED.RData")

UMAP plot for Hannifa healthy skin data

DimPlot(healthy_skin.sc_data,raster = FALSE)

pdf(file = "UMAP_ONLY_HEALTHY_SAMPLES_(HANNIFA_PAPER).pdf",width = 16,height=8)
print(DimPlot(healthy_skin.sc_data_subset,pt.size = 1.2,shuffle = TRUE,raster=FALSE))
dev.off()

Identify marker genes for the pre-assigned clustering results with cell type labels

Idents(healthy_skin.sc_data) <- "final_clustering"
hnf_skin_scRNA.markers <- FindAllMarkers(healthy_skin.sc_data,assay = "RNA",logfc.threshold = 0.25,max.cells.per.ident=400,only.pos = TRUE)

Load MARKER GENES from the saved file for later use

hnf_skin_scRNA.markers <- read.csv("../../scRNA_DATA/HANNIFA_MARKERS.csv") %>% filter(p_val_adj<=0.05) %>% group_by(cluster) %>% top_n(n =300,wt = avg_log2FC) %>% filter(avg_log2FC>0.25)

#SC-DATA - cell counts

hnf_data_counts <- as.data.frame(table(healthy_skin.sc_data@meta.data$final_clustering))
write.csv(hnf_data_counts,file="HEALTHY_SAMPLES_HANNIFA_DATA_COUNTS.csv")

Intersecting genes between ST and scRNA data.

##INTERSECT GENES BETWEEN scRNA and Spatial data
st.genes <- unique(rownames(new.skin.combined@assays$Spatial@counts))
sc.genes <- unique(rownames(healthy_skin.sc_data@assays$RNA@counts))
all.genes.scrna_and_spt <- unique(intersect(sc.genes,st.genes))

MIA FOR HANNIFA PAPER - healthy samples data only

MIA_results_PART_2 <- MIA(total_genes = length(all.genes.scrna_and_spt),single_cell.markers = hnf_skin_scRNA.markers,spatial.markers = filtered_spatial_markers)

E.data_PART_2 <- MIA_results_PART_2 %>% column_to_rownames("cluster")
E.data_PART_2 <- E.data_PART_2[,order(colnames(E.data_PART_2))]
is.na(E.data_PART_2) <- do.call(cbind,lapply(E.data_PART_2, is.infinite))
pheatmap(E.data_PART_2,cluster_cols = FALSE,cluster_rows = FALSE,fontsize=15,color = col.pal)

For Immune populations

## IMMUNE CELLS
immune_only.E.data_PART_2 <- E.data_PART_2[,c("DC1","DC2","Macro_1","Macro_2","moDC_1","moDC_2","moDC_3","MigDC","Mono","Inf_mono","ILC1_NK","ILC1_3","ILC2","Tc","Th","Treg","Mast_cell","Plasma","NK")]
pdf(file = "MIA_regions_(HANNIFA_PAPER)_IMMUNE_CELLS.pdf",width = 10,height = 5)
print(pheatmap(immune_only.E.data_PART_2,cluster_cols = FALSE,cluster_rows = FALSE,fontsize=15,color = col.pal))
dev.off()

For non-immune populations

non_immune_only.E.data_PART_2 <- E.data_PART_2 %>% dplyr::select(-c("DC1","DC2","Macro_1","Macro_2","moDC_1","moDC_2","moDC_3","MigDC","Mono","Inf_mono","ILC1_NK","ILC1_3","ILC2","Tc","Th","Treg","Mast_cell","Plasma","NK","nan"))
pdf(file = "MIA_regions_(HANNIFA_PAPER)_NON_IMMUNE_CELLS.pdf",width = 10,height = 5)
print(pheatmap(non_immune_only.E.data_PART_2,cluster_cols = FALSE,cluster_rows = FALSE,fontsize=15,color = col.pal))
dev.off()

#VIOLIN PLOT FOR UMIs and SEURAT CLUSTERS
pdf(file = "VIOLIN_PLOT_UMIs_and_CLUSTERS_HEALTHY_SAMPLES.pdf",width = 16,height=12)
print(VlnPlot(new.skin.combined,group.by = "seurat_clusters",features = "nCount_Spatial",cols = custom_colors))
dev.off()

# LOG TRANSFORMED UMI COUNTS
pdf(file = "VIOLIN_PLOT_UMIs_and_CLUSTERS_HEALTHY_SAMPLES(LOG_TRANSFORMED).pdf",width = 16,height=12)
print(VlnPlot(new.skin.combined,group.by = "seurat_clusters",features = "nCount_Spatial",cols = custom_colors,log = TRUE))
dev.off()

# SPOT FREQUENCY TABLE - BY SPATIAL REGION
spatial.freq <- as.data.frame(table(new.skin.combined@meta.data$Spatial.regions))
write.csv(spatial.freq,file="SPATIAL_REGIONS_CLUSTERS_SPOT_COUNTS.csv")

QC FIGURES - Replotting

Supplementary figure - 2

# HEALTHY FEMALE SKIN 2
#REMOVE SPOTS
HEALTHY.Female1.s3 <- Load10X_Spatial(data.dir ="../../SEVENTH RUN/ST-HF-1/",slice="HealthyFemale-1-S3")

## HEALTHY FEMALE
HEALTHY.Female2.s1 <- Load10X_Spatial(data.dir ="../../SEVENTH RUN/ST-HF-2E/",slice="HealthyFemale-2-S1")
HEALTHY.Female2.s2 <- Load10X_Spatial(data.dir ="../../SEVENTH RUN/ST-HF-2F/",slice="HealthyFemale-2-S2")

#HEALATHY MALE
HEALTHY.Male.s2 <- Load10X_Spatial(data.dir ="../../SEVENTH RUN/ST-HM-1/",slice="HealthyMale-1-S2")

remove.spots <- read.csv(file="../../CLOUPE_FILES/SEVENTH_RUN/HF-1-SP/remove.csv")
subset_spots <- Cells(HEALTHY.Female1.s3)[which((!(rownames(HEALTHY.Female1.s3@meta.data) %in% remove.spots$Barcode)))]
HEALTHY.Female1.s3.clean <- subset(HEALTHY.Female1.s3,cells=subset_spots)

pdf(width = 8,height=10,file = "HV2_SP_SPATIAL_PLOT(UNFILTERED).pdf")
print(SpatialDimPlot(HEALTHY.Female1.s3.clean,pt.size.factor = 2,group.by = "orig.ident"))
dev.off()

pdf(width = 8,height=10,file = "HV2_SP_SPATIAL_PLOT(FILTERED).pdf")
print(SpatialDimPlot(new.skin.combined,images = "HealthyFemale.1.S3",pt.size.factor = 2,group.by = "orig.ident"))
dev.off()

remove.spots <- read.csv(file="../../CLOUPE_FILES/SEVENTH_RUN/HF-2-SE/Remove.csv")
subset_spots <- Cells(HEALTHY.Female2.s1)[which((!(rownames(HEALTHY.Female2.s1@meta.data) %in% remove.spots$Barcode)))]
HEALTHY.Female2.s1.clean <- subset(HEALTHY.Female2.s1,cells=subset_spots)

pdf(width = 8,height=10,file = "HV3_SE_SPATIAL_PLOT(UNFILTERED).pdf")
print(SpatialDimPlot(HEALTHY.Female2.s1.clean,pt.size.factor = 2,group.by = "orig.ident"))
dev.off()

pdf(width = 8,height=10,file = "HV3_SE_SPATIAL_PLOT(FILTERED).pdf")
print(SpatialDimPlot(new.skin.combined,images = "HealthyFemale.2.S1",pt.size.factor = 2,group.by = "orig.ident"))
dev.off()

remove.spots <- read.csv(file="../../CLOUPE_FILES/SEVENTH_RUN/HF-2-SP/Remove.csv")
subset_spots <- Cells(HEALTHY.Female2.s2)[which((!(rownames(HEALTHY.Female2.s2@meta.data) %in% remove.spots$Barcode)))]
HEALTHY.Female2.s2.clean <- subset(HEALTHY.Female2.s2,cells=subset_spots)

pdf(width = 8,height=10,file = "HV3_SP_SPATIAL_PLOT(UNFILTERED).pdf")
print(SpatialDimPlot(HEALTHY.Female2.s2.clean,pt.size.factor = 2,group.by = "orig.ident"))
dev.off()

pdf(width = 8,height=10,file = "HV3_SP_SPATIAL_PLOT(FILTERED).pdf")
print(SpatialDimPlot(new.skin.combined,images = "HealthyFemale.2.S2",pt.size.factor = 2,group.by = "orig.ident"))
dev.off()

remove.spots <- read.csv(file="../../CLOUPE_FILES/SEVENTH_RUN/HM-SP-2/remove.csv")
subset_spots <- Cells(HEALTHY.Male.s2)[which((!(rownames(HEALTHY.Male.s2@meta.data) %in% remove.spots$Barcode)))]
HEALTHY.Male.s2.clean <- subset(HEALTHY.Male.s2,cells=subset_spots)

pdf(width = 8,height=10,file = "HV1_SP_R2_SPATIAL_PLOT(UNFILTERED).pdf")
print(SpatialDimPlot(HEALTHY.Male.s2.clean,pt.size.factor = 2.2,group.by = "orig.ident"))
dev.off()

pdf(width = 8,height=10,file = "HV1_SP_R2_SPATIAL_PLOT(FILTERED).pdf")
print(SpatialDimPlot(new.skin.combined,images = "HealthyMale.1.S2",pt.size.factor = 2.2,group.by = "orig.ident"))
dev.off()

Import the spots to be removed from loupe browser.

remove.spots <- read.csv(file="../../FIFTH_RUN/SAMPLE D1/REMOVE.csv")
subset_spots <- Cells(HEALTHY.Female1.s2)[which((!(rownames(HEALTHY.Female1.s2@meta.data) %in% remove.spots$Barcode)))]
HEALTHY.Female1.s2.clean <- subset(HEALTHY.Female1.s2,cells=subset_spots)

SpatialDimPlot(HEALTHY.Female.s2.clean)
SpatialDimPlot(HEALTHY.Female1.s2.clean)

pdf(width = 8,height=10,file = "HV2_SP_R2_SPATIAL_PLOT(UNFILTERED).pdf")
print(SpatialDimPlot(HEALTHY.Female1.s2.clean,pt.size.factor = 2.2,group.by = "orig.ident"))
dev.off()

pdf(width = 8,height=10,file = "HV2_SP_R2_SPATIAL_PLOT(FILTERED).pdf")
print(SpatialDimPlot(new.skin.combined,images = "HealthyFemale.1.S2",pt.size.factor = 2.2,group.by = "orig.ident"))
dev.off()

AVERAGE GENE COUNTS

Idents(new.skin.combined) <- "seurat_clusters"

clusters <- unique(new.skin.combined@meta.data$seurat_clusters) %>% as.vector()

for(x in clusters){
  avg_exp.cluster <- AverageExpression(new.skin.combined,slot = "counts",assay="Spatial") %>% as.data.frame() %>% dplyr::select(paste("Spatial.",x,sep = "")) %>% filter(x>0) %>% write.csv(file = paste("CLUSTER_",x,"_GENE_COUNTS.csv",sep = ""))
}

avg_exp.all_clusters <- AverageExpression(new.skin.combined,slot = "counts",assay="Spatial") %>% as.data.frame()
write.csv(avg_exp.all_clusters,file="ALL_CLUSTERS_GENE_COUNTS.csv")

avg_exp.cluster_4 <- AverageExpression(new.skin.combined,slot = "counts",assay="Spatial") %>% as.data.frame() %>% dplyr::select("Spatial.4.Hypodermis") %>% filter(Spatial.4.Hypodermis>0)
write.csv(avg_exp.cluster_4,file="CLUSTER_4_GENE_COUNTS.csv")
avg_exp.cluster_2 <- AverageExpression(new.skin.combined,slot = "counts",assay="Spatial") %>% as.data.frame() %>% dplyr::select("Spatial.2.Dermis.adipose") %>% filter(Spatial.2.Dermis.adipose>0)
write.csv(avg_exp.cluster_2,file="CLUSTER_2_GENE_COUNTS.csv")

HANNIFA DATA - NEW UMAP

healthy_skin.sc_data <- FindNeighbors(healthy_skin.sc_data, dims = 1:40)
healthy_skin.sc_data <- FindClusters(healthy_skin.sc_data, verbose = FALSE,resolution = 1)
healthy_skin.sc_data <- RunUMAP(healthy_skin.sc_data, dims = 1:40)

Idents(healthy_skin.sc_data) <- "final_clustering"
healthy_skin.sc_data <- subset(healthy_skin.sc_data,idents = c("nan"),invert=TRUE)

pdf(width = 12,height=8,file = "UMAP_HANNIFA_DATA_HEALTHY_SAMPLES_ONLY.pdf")
DimPlot(healthy_skin.sc_data,pt.size = 3.5,raster=FALSE)
dev.off()

SINGLE CELL DATA - ENRICHMENT PLOTS (USING ONLY ONE HEALTHY SAMPLE)

HEALTHY.Female1.s2.clean<- st_filter_by_genes(st.data = HEALTHY.Female1.s2.clean,x = 200)
HEALTHY.Female1.s2.clean <- SCTransform(HEALTHY.Female1.s2.clean, verbose = FALSE,assay="Spatial")
HEALTHY.Female1.s2.clean <- RunPCA(HEALTHY.Female1.s2.clean)

#travis_NM_skin_scRNA <- readRDS(file="../../TRAVIS_scRNA_DATA/travis_NM_skin_scRNA.RDS")
# RUN THIS ON BIG PURPLE 
skin_reference.3 <- SCTransform(travis_NM_skin_scRNA, ncells = 3000, verbose = FALSE) %>% RunPCA(verbose = FALSE) %>% 
    RunUMAP(dims = 1:40)

iteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachedNaNs producediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachedNaNs producediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachedNaNs producediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachedNaNs producediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachedNaNs producediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachedNaNs producediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachedNaNs producediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachedNaNs producediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachedNaNs producediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachedNaNs producediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachedNaNs producediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachedNaNs producediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachedNaNs producediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachediteration limit reachedThe default method for RunUMAP has changed from calling Python UMAP via reticulate to the R-native UWOT using the cosine metric
To use Python UMAP via reticulate, set umap.method to 'umap-learn' and metric to 'correlation'
This message will be shown once per session15:43:15 UMAP embedding parameters a = 0.9922 b = 1.112
15:43:15 Read 4840 rows and found 40 numeric columns
15:43:15 Using Annoy for neighbor search, n_neighbors = 30
15:43:15 Building Annoy index with metric = cosine, n_trees = 50
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
15:43:15 Writing NN index file to temp file /var/folders/8p/2clwbcfn719f5l1kjsyx4v6s2r9kgj/T//Rtmpovk9y0/filefd7577fa79e5
15:43:15 Searching Annoy index using 1 thread, search_k = 3000
15:43:16 Annoy recall = 100%
15:43:17 Commencing smooth kNN distance calibration using 1 thread
15:43:18 Initializing from normalized Laplacian + noise
15:43:19 Commencing optimization for 500 epochs, with 191228 positive edges
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
15:43:25 Optimization finished

anchors <- FindTransferAnchors(reference = skin_reference.3, query = HEALTHY.Female1.s2.clean, normalization.method = "SCT")

Normalizing query using reference SCT model
Adding image data that isn't associated with any assay presentPerforming PCA on the provided reference using 2702 features as input.
Projecting cell embeddings
Finding neighborhoods
Finding anchors
    Found 530 anchors

predictions.assay <- TransferData(anchorset = anchors, refdata = skin_reference.3$Specific_CellType, prediction.assay = TRUE, 
    weight.reduction = HEALTHY.Female1.s2.clean[["pca"]],dims = 1:40)

Finding integration vectors
Finding integration vector weights
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Predicting cell labels

HEALTHY.Female1.s2.clean[["predictions_travis_data"]] <- predictions.assay

library(RColorBrewer)
SpatialColors <- colorRampPalette(colors = rev(x = brewer.pal(n = 11, name = "Spectral")))
cols <- SpatialColors(n = 100)

SpatialFeaturePlot(HEALTHY.Female1.s2.clean,features = c("Melanocyte"), pt.size.factor = 2.5)+ scale_fill_gradientn(colors=cols,limits = c(0,1.000001))

Could not find Melanocyte in the default search locations, found in predictions_travis_data assay instead`guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.Scale for 'fill' is already present. Adding another scale for 'fill', which will replace the existing scale.

SpatialFeaturePlot(HEALTHY.Female1.s2.clean,features = c("VSMC-1"), pt.size.factor = 2.5)+ scale_fill_gradientn(colors=cols,limits = c(0,1.000001))

Could not find VSMC-1 in the default search locations, found in predictions_travis_data assay instead`guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.Scale for 'fill' is already present. Adding another scale for 'fill', which will replace the existing scale.

SpatialFeaturePlot(HEALTHY.Female1.s2.clean,features = c("HairFollicle"), pt.size.factor = 2.5)+ scale_fill_gradientn(colors=cols,limits = c(0,1.000001))

Could not find HairFollicle in the default search locations, found in predictions_travis_data assay instead`guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.Scale for 'fill' is already present. Adding another scale for 'fill', which will replace the existing scale.

SpatialFeaturePlot(HEALTHY.Female1.s2.clean,features = c("Keratinocyte-5"), pt.size.factor = 2.5)+ scale_fill_gradientn(colors=cols,limits = c(0,1.000001))

Could not find Keratinocyte-5 in the default search locations, found in predictions_travis_data assay instead`guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.Scale for 'fill' is already present. Adding another scale for 'fill', which will replace the existing scale.

SpatialFeaturePlot(HEALTHY.Female1.s2.clean,features = c("Sebocyte"), pt.size.factor = 2.5)+ scale_fill_gradientn(colors=cols,limits = c(0,1.000001))

Could not find Sebocyte in the default search locations, found in predictions_travis_data assay instead`guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.Scale for 'fill' is already present. Adding another scale for 'fill', which will replace the existing scale.